Factor analysis
library(readxl)
mydata <- read_xlsx("./Purchasing decision.xlsx")
mydata <- as.data.frame(mydata)
head(mydata)## ID X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12
## 1 1 5 5 1 3 4 5 1 3 5 3 3 1
## 2 2 3 2 1 2 1 2 2 2 3 2 3 1
## 3 3 5 4 1 4 4 4 1 3 5 3 4 2
## 4 4 2 3 1 2 3 3 1 3 2 3 2 3
## 5 5 4 4 2 3 5 3 3 5 4 5 2 5
## 6 6 4 3 3 2 2 4 1 2 5 2 4 1All variables are measured on Likert scale (1: Not important at all, 5: Extremely important)
- ID
- X1: Reliable
- X2: Attractive
- X3: Practical
- X4: Strong
- X5: Exciting
- X6: Innovative
- X7: Robust
- X8: Fun
- X9: Build quality
- X10:Powerful
- X11:Versatile
- X12:Sporty
mydata_FA <- mydata[c(2:13)]
desc_stats <- data.frame(
N = apply(mydata_FA, 2, function(x) sum(!is.na(x))),
Missing = apply(mydata_FA, 2, function(x) sum(is.na(x))),
Mean = apply(mydata_FA, 2, mean, na.rm = TRUE),
Variance = apply(mydata_FA, 2, var, na.rm = TRUE),
Min = apply(mydata_FA, 2, min, na.rm = TRUE),
Max = apply(mydata_FA, 2, max, na.rm = TRUE)
)
round(desc_stats, 2)## N Missing Mean Variance Min Max
## X1 303 0 3.67 1.04 1 5
## X2 303 0 3.49 1.21 1 5
## X3 303 0 2.66 1.95 1 5
## X4 302 1 3.60 1.01 1 5
## X5 303 0 3.44 1.50 1 5
## X6 303 0 3.70 0.71 1 5
## X7 301 2 2.52 1.76 1 5
## X8 302 1 3.48 1.57 1 5
## X9 302 1 3.75 0.95 1 5
## X10 302 1 3.44 1.37 1 5
## X11 302 1 3.09 1.17 1 5
## X12 299 4 3.03 1.72 1 5mydata <- mydata[complete.cases(mydata), ]
mydata_FA <- mydata[c(2:13)]
nrow(mydata_FA)## [1] 296R <- cor(mydata_FA) #Saving correlation matrixlibrary(GGally)
ggpairs(
mydata_FA,
upper = list(continuous = wrap("points", size = 1)),
lower = list(continuous = wrap("cor", size = 2.5), colour = "black"),
diag = list(continuous = wrap("densityDiag", linewidth = 0.6))) +
theme(
axis.text.x = element_text(size = 5),
axis.text.y = element_text(size = 5),
strip.text.x = element_text(size = 6),
strip.text.y = element_text(size = 5),
axis.title = element_blank())
# Compute clustering order from correlation matrix
hc <- hclust(as.dist(1 - abs(R)), method = "ward.D2")
# Reorder correlation matrix
R_ord <- R[hc$order, hc$order]mydata_FA_ord <- mydata_FA[, hc$order]
ggpairs(
mydata_FA_ord,
upper = list(continuous = wrap("points", size = 1)),
lower = list(continuous = wrap("cor", size = 2.5), colour = "black"),
diag = list(continuous = wrap("densityDiag", linewidth = 0.6))) +
theme(
axis.text.x = element_text(size = 5),
axis.text.y = element_text(size = 5),
strip.text.x = element_text(size = 6),
strip.text.y = element_text(size = 5),
axis.title = element_blank())
det(R) #Determinant of correlation matrix## [1] 0.00095629library(psych)
cortest.bartlett(R, n = nrow(mydata_FA)) #Bartlett's test of sphericity## $chisq
## [1] 2017.369
##
## $p.value
## [1] 0
##
## $df
## [1] 66library(psych)
KMO(R) #KMO and MSA statistics## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = R)
## Overall MSA = 0.81
## MSA for each item =
## X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12
## 0.82 0.85 0.69 0.89 0.86 0.76 0.67 0.88 0.70 0.86 0.84 0.82library(psych)
fa.parallel(
mydata_FA,
fa = "fa",
fm = "minres",
cor = "cor",
n.iter = 1000,
sim = TRUE,
error.bars = TRUE
)## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Parallel analysis suggests that the number of factors = 3 and the number of components = NAlibrary(psych)
library(GPArotation)
factors <- fa_ordered(
mydata_FA,
nfactors = 3,
covar = FALSE,
fm = "minres",
rotate = "oblimin",
scores = "regression")
print.psych(
factors,
cut = 0.40,
sort = TRUE)## Factor Analysis using method = minres
## Call: psych::fa(r = r, nfactors = nfactors, rotate = "oblimin", scores = "regression",
## covar = FALSE, fm = "minres")
## Standardized loadings (pattern matrix) based upon correlation matrix
## item F1 F2 F3 h2 u2 com
## X5 5 0.87 0.77 0.23 1.0
## X8 8 0.86 0.75 0.25 1.0
## X10 10 0.78 0.65 0.35 1.2
## X12 12 0.72 0.71 0.29 1.8
## X4 4 0.70 0.61 0.39 1.3
## X2 2 0.61 0.41 0.66 0.34 2.2
## X7 7 0.90 0.82 0.18 1.0
## X3 3 0.88 0.78 0.22 1.0
## X11 11 0.55 0.41 0.59 2.0
## X9 9 0.82 0.66 0.34 1.1
## X6 6 0.74 0.58 0.42 1.1
## X1 1 0.64 0.43 0.57 1.1
##
## F1 F2 F3
## SS loadings 3.59 2.19 2.04
## Proportion Var 0.30 0.18 0.17
## Cumulative Var 0.30 0.48 0.65
## Proportion Explained 0.46 0.28 0.26
## Cumulative Proportion 0.46 0.74 1.00
##
## With factor correlations of
## F1 F2 F3
## F1 1.00 0.07 0.09
## F2 0.07 1.00 -0.07
## F3 0.09 -0.07 1.00
##
## Mean item complexity = 1.3
## Test of the hypothesis that 3 factors are sufficient.
##
## df null model = 66 with the objective function = 6.95 with Chi Square = 2017.37
## df of the model are 33 and the objective function was 0.28
##
## The root mean square of the residuals (RMSR) is 0.02
## The df corrected root mean square of the residuals is 0.03
##
## The harmonic n.obs is 296 with the empirical chi square 10.02 with prob < 1
## The total n.obs was 296 with Likelihood Chi Square = 80.8 with prob < 6.9e-06
##
## Tucker Lewis Index of factoring reliability = 0.951
## RMSEA index = 0.07 and the 90 % confidence intervals are 0.051 0.09
## BIC = -106.98
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## F1 F2 F3
## Correlation of (regression) scores with factors 0.96 0.95 0.91
## Multiple R square of scores with factors 0.92 0.90 0.83
## Minimum correlation of possible factor scores 0.84 0.80 0.67indicator_labels <- c(
X1 = "Reliable", X2 = "Attractive", X3 = "Practical", X4 = "Strong",
X5 = "Exciting", X6 = "Innovative", X7 = "Robust", X8 = "Fun",
X9 = "High build quality", X10 = "Powerful", X11 = "Versatile", X12 = "Sporty")
factor_labels <- c(
F1 = "Performance",
F2 = "Functionality",
F3 = "Quality and Reliability")
loadings_diagram <- as.matrix(factors$loadings)
rownames(loadings_diagram) <- indicator_labels[rownames(loadings_diagram)]
colnames(loadings_diagram) <- factor_labels[colnames(loadings_diagram)]
Phi_diagram <- factors$Phi
dimnames(Phi_diagram) <- list(colnames(loadings_diagram), colnames(loadings_diagram))
psych::fa.diagram(
loadings_diagram,
Phi = Phi_diagram,
cut = 0.40,
simple = TRUE,
main = "Factor Analysis")
R_obs <- factors$r
round(R_obs, 2)## X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12
## X1 1.00 0.40 -0.06 0.16 0.17 0.45 0.00 0.09 0.55 0.23 0.13 -0.09
## X2 0.40 1.00 -0.19 0.44 0.64 0.47 -0.16 0.47 0.26 0.55 -0.13 0.21
## X3 -0.06 -0.19 1.00 0.32 0.00 0.02 0.80 0.15 0.02 -0.13 0.44 0.34
## X4 0.16 0.44 0.32 1.00 0.59 0.22 0.30 0.63 0.09 0.58 0.07 0.55
## X5 0.17 0.64 0.00 0.59 1.00 0.23 0.04 0.74 -0.02 0.68 -0.14 0.57
## X6 0.45 0.47 0.02 0.22 0.23 1.00 0.05 0.13 0.58 0.19 0.25 -0.12
## X7 0.00 -0.16 0.80 0.30 0.04 0.05 1.00 0.15 0.09 -0.13 0.46 0.38
## X8 0.09 0.47 0.15 0.63 0.74 0.13 0.15 1.00 -0.08 0.62 -0.11 0.68
## X9 0.55 0.26 0.02 0.09 -0.02 0.58 0.09 -0.08 1.00 0.10 0.30 -0.22
## X10 0.23 0.55 -0.13 0.58 0.68 0.19 -0.13 0.62 0.10 1.00 -0.20 0.47
## X11 0.13 -0.13 0.44 0.07 -0.14 0.25 0.46 -0.11 0.30 -0.20 1.00 0.00
## X12 -0.09 0.21 0.34 0.55 0.57 -0.12 0.38 0.68 -0.22 0.47 0.00 1.00R_rep <- factors$model
round(R_rep, 2)## X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12
## X1 0.43 0.39 -0.06 0.17 0.19 0.50 -0.02 0.09 0.50 0.21 0.14 -0.10
## X2 0.39 0.66 -0.20 0.44 0.60 0.43 -0.17 0.49 0.30 0.59 -0.11 0.23
## X3 -0.06 -0.20 0.78 0.29 0.01 0.01 0.79 0.15 0.04 -0.12 0.45 0.36
## X4 0.17 0.44 0.29 0.61 0.62 0.22 0.31 0.64 0.09 0.52 0.07 0.57
## X5 0.19 0.60 0.01 0.62 0.77 0.21 0.02 0.74 0.02 0.69 -0.15 0.57
## X6 0.50 0.43 0.01 0.22 0.21 0.58 0.06 0.11 0.59 0.23 0.21 -0.09
## X7 -0.02 -0.17 0.79 0.31 0.02 0.06 0.82 0.16 0.09 -0.10 0.47 0.36
## X8 0.09 0.49 0.15 0.64 0.74 0.11 0.16 0.75 -0.08 0.64 -0.10 0.66
## X9 0.50 0.30 0.04 0.09 0.02 0.59 0.09 -0.08 0.66 0.05 0.30 -0.25
## X10 0.21 0.59 -0.12 0.52 0.69 0.23 -0.10 0.64 0.05 0.65 -0.19 0.44
## X11 0.14 -0.11 0.45 0.07 -0.15 0.21 0.47 -0.10 0.30 -0.19 0.41 -0.02
## X12 -0.10 0.23 0.36 0.57 0.57 -0.09 0.36 0.66 -0.25 0.44 -0.02 0.71R_res <- factors$residual
round(R_res,2)## X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12
## X1 0.57 0.01 0.00 -0.01 -0.02 -0.04 0.02 0.00 0.04 0.02 -0.02 0.01
## X2 0.01 0.34 0.01 0.00 0.04 0.04 0.01 -0.01 -0.04 -0.04 -0.02 -0.02
## X3 0.00 0.01 0.22 0.03 0.00 0.01 0.01 0.00 -0.02 -0.01 -0.01 -0.03
## X4 -0.01 0.00 0.03 0.39 -0.03 -0.01 -0.01 -0.01 0.00 0.06 0.00 -0.02
## X5 -0.02 0.04 0.00 -0.03 0.23 0.03 0.01 0.01 -0.03 -0.02 0.01 0.00
## X6 -0.04 0.04 0.01 -0.01 0.03 0.42 -0.01 0.02 -0.01 -0.03 0.03 -0.03
## X7 0.02 0.01 0.01 -0.01 0.01 -0.01 0.18 -0.01 0.00 -0.02 -0.01 0.02
## X8 0.00 -0.01 0.00 -0.01 0.01 0.02 -0.01 0.25 0.00 -0.02 0.00 0.02
## X9 0.04 -0.04 -0.02 0.00 -0.03 -0.01 0.00 0.00 0.34 0.05 0.00 0.02
## X10 0.02 -0.04 -0.01 0.06 -0.02 -0.03 -0.02 -0.02 0.05 0.35 -0.02 0.03
## X11 -0.02 -0.02 -0.01 0.00 0.01 0.03 -0.01 0.00 0.00 -0.02 0.59 0.02
## X12 0.01 -0.02 -0.03 -0.02 0.00 -0.03 0.02 0.02 0.02 0.03 0.02 0.29Residuals
residual_values <- factors$residual[upper.tri(factors$residual)]
high_residuals <- abs(residual_values) > 0.05
mean(high_residuals)## [1] 0.01515152Saving factor values
scores_ordered <- factors$scores # ze F1, F2, F3 v pravilnem vrstnem redu
mydata$F1_Performance <- scores_ordered[, "F1"]
mydata$F2_Functionality <- scores_ordered[, "F2"]
mydata$F3_Quality_and_Reliability <- scores_ordered[, "F3"]
head(mydata[, c("ID", "F1_Performance", "F2_Functionality",
"F3_Quality_and_Reliability")], 3)## ID F1_Performance F2_Functionality F3_Quality_and_Reliability
## 1 1 -0.2226898 -1.3013388 1.729783
## 2 2 -1.7660244 -0.6734816 -1.038851
## 3 3 -0.1053196 -1.0207023 1.198150head(mydata[1, c(-1)])## X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12 F1_Performance F2_Functionality
## 1 5 5 1 3 4 5 1 3 5 3 3 1 -0.2226898 -1.301339
## F3_Quality_and_Reliability
## 1 1.729783Cronbach alfa
Performance <- mydata[, c("X5", "X8", "X10", "X12", "X4","X2")] #Select only indicators that measure factor Performance
library(psych)
alpha(Performance, #Data
check.keys = FALSE) ##
## Reliability analysis
## Call: alpha(x = Performance, check.keys = FALSE)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.88 0.88 0.89 0.56 7.7 0.01 3.4 0.95 0.58
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.86 0.88 0.9
## Duhachek 0.86 0.88 0.9
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## X5 0.84 0.84 0.84 0.52 5.4 0.014 0.0181 0.55
## X8 0.84 0.85 0.85 0.53 5.6 0.014 0.0175 0.56
## X10 0.86 0.86 0.86 0.55 6.2 0.013 0.0230 0.58
## X12 0.88 0.88 0.87 0.59 7.3 0.011 0.0082 0.61
## X4 0.87 0.87 0.87 0.56 6.5 0.012 0.0235 0.59
## X2 0.89 0.89 0.87 0.61 7.9 0.010 0.0061 0.61
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## X5 296 0.89 0.88 0.87 0.82 3.4 1.2
## X8 296 0.88 0.87 0.85 0.80 3.5 1.3
## X10 296 0.81 0.82 0.76 0.72 3.4 1.2
## X12 296 0.74 0.73 0.67 0.61 3.0 1.3
## X4 296 0.78 0.79 0.73 0.69 3.6 1.0
## X2 296 0.68 0.69 0.62 0.55 3.5 1.1
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## X5 0.09 0.14 0.23 0.32 0.22 0
## X8 0.09 0.13 0.22 0.30 0.25 0
## X10 0.09 0.11 0.27 0.34 0.19 0
## X12 0.16 0.21 0.26 0.21 0.17 0
## X4 0.03 0.11 0.29 0.38 0.19 0
## X2 0.05 0.13 0.28 0.34 0.19 0